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' A recently introduced particle-based model for fluid dynamics with effective excluded volume interactions is analyzed 
in detail. The interactions are modeled by means of stochastic multiparticle collisions which are biased and depend 
on local velocities and densities. Momentum and energy are exactly conserved locally. The isotropy and relaxation 
to equilibrium are analyzed and measured. It is shown how a discrete-time projection operator technique can be 
used to obtain Green-Kubo relations for the transport coefficients. Because of a large viscosity no long-time tails in 
the velocity auto-correlation and stress correlation functions were seen. Strongly reduced self-diffusion due to caging 
and an order/disorder transition is found at high collision frequency, where clouds consisting of at least four particles 
form a cubic phase. These structures were analyzed by measuring the pair-correlation function above and below the 
transition. Finally, the algorithm is extended to binary mixtures which phase-separate above a critical collision rate, 
g ; PACS number(s): 47.11. +j, 05.40. +j, 02.70.Ns 
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I. INTRODUCTION 



The efficient modeling of the hydrodynamics of complex liquids such as colloidal suspensions and microemulsions 
is a challenging task due to the interplay between flow and thermodynamic properties and the importance of thermal 
noise. The desire to simplify the description of the solvent degrees of freedom led to the development of mesoscale 
simulation methods, which are essentially novel ways of solving the equations of fluctuating hydrodynamics. 

In particular, a recently introduced simulation technique of this type — often called Stochastic Rotation Dynamics 
method (SRD) 1] or Multi-Particle Collision Dynamics — is a very promising tool for such a coarse-grained 
description of a fluctuating fluid 0, . The method is based on so-called fluid particles with continuous positions and 
velocities, which follow simple dynamic rules of streaming and collision. In this sense, SRD is closer to the microscopic 
reality than conventional Navier-Stokes solvers. On the one hand, this gives rise to unlimited numerical stability, on 
the other hand there are certain restrictions on the transport coefficients which depend on the details of the dynamics. 
For instance, the viscosity cannot be chosen to be arbitrarily small or large at reasonable computational efficiency 
Due to the simplicity of the SRD-rules, many analytical calculations could be performed |5l|6| which are very hard to 
do for other particle-based methods. So far, SRD has been used to study sedimentation [jUll, colloids and polymers 
in flow [U ^3 ' vesicle dynamics ^lj , reacting flow 0] > an d microemulsions . 

The fluid particles of the original SRD method have an ideal gas equation of state. Hence, they are very compressible 
and the speed of sound, c s , is low. In order to have negliglible compressibility effects as in real liquids the Mach number 
has to be kept small, which means that there are limits for the flow speed in the simulation. Recently, we showed how 
SRD can be modified to achieve a larger speed of sound and a larger viscosity 0, which makes simulations of dense 
gases and liquids more efficient and realistic. This was done in a thermodynamically consistent way by introducing 
generalized excluded volume interactions among the fluid particles. That means, the kinetic energy is still locally 
conserved and no complications such as temperature drifts are encountered. This was a first step towards a consistent 
SRD-model for more general liquids with additional attractive interactions and the possibility of a liquid-gas phase 
transition. 

While (to our knowledge) this has been the first SRD model with a non-ideal equation of state, there have been 
several attempts to generalize other particle methods to non-ideal fluids. The first one is described by Alexander et al, 
[l7j . and is based on a modification of the original Direct Simulation Monte Carlo (DSMC) or Bird's algorithm [Isj . 
It is called the Consistent Boltzmann algorithm (CBA) |17ill9j . The authors labeled the original DSMC algorithm as 
"inconsistent" in the sense that it has the transport properties of a hard-sphere gas, but obeys an ideal gas equation of 
state; CBA has the correct hard-sphere equation of state. However, back-scattering events connected with structural 
effects are absent in the model leading to too high self-diffusion at higher densities. Furthermore, the sound speed 
determined from the equal-time density fluctuations are not in agreement with the value obtained directly from the 
equation of state. As seen from Fig. 2 and 3 in the analytical expressions from Enskog-theory for the transport 
coefficients are significantly lower than the simulation results. 

Another particle method very close to DSMC explicitly aimed to solve the Enskog equation is given in Rcfs. 
jiol Elf . The method correctly reproduced the transport properties of the Enskog gas, but has the unpleasant feature 
of conserving momentum and energy only in a statistical sense and not in a single collision. To overcome these 
deficiencies a modification of the DSMC method was proposed to solve the Enskog-equation by Frezzotti [i^] ■ The 
main idea is the same as in 21] : due to the finite size of a particle, binary collisons do not just take place between 
particles of the same cell, but also particles from adjacent cells can participate, but now energy and momentum are 
conserved in every collision. 

The original SRD method has a similiar issue as DSMC: The transport properties can be tuned to that of an 
interacting gas, while the equation of state remains ideal. Hence, a generalization which leads to a tunable equation 
of state is highly desirable. This also allows us to observe phase transitions and critical phenomena, and the method 
can be easily extended to multi-component systems. Furthermore, it was shown that the generalized SRD is "more 
consistent" than CBA since independent measurements of the density fluctuations, the pressure and the speed of 
sound are consistent with each other and agree with theory • 

In CBA, excluded volume effects are modeled by introducing an additional displacement of the particles away from 
each other after every binary collision. This displacement is in addition to the free streaming displacement, TVi, 
given by the time step r and the actual velocity Vi of particle i. The idea behind this is that particles of finite size 
collide earlier than point-like particles, and hence would be further apart after time r. It is likely that this enhanced 
particle relocation leads to an additional diffusive term in the continuity equation for the density, which could be 
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the reason for thermodynamically inconsistent density fluctuations. As discussed in [15( one has to be careful with 
respect to thermodynamic consistency; often it is not obvious why a certain algorithm fails to produce the correct 
thermal fluctuations. Therefore, in generalizing SRD we followed a different strategy than Alexander et al. [l7|: Our 
main idea is to coarse-grain the collisions between hard-spheres. This leads to additional effects such as caging and 
crystallization, and thermodynamically consistent fluctuations can be obtained. 

Instead of costly binary collisions, which is the core of any DSMC-code, imagine two clouds of particles containing 
Mi and Mi particles, and having center of mass velocities u\, 1*2 respectively. The volume of the clouds is assumed to 
be equal and constant. On average, the momentum transfer between these clouds will be high if they are very dense 
and move at high speed into each other. Hence, the average interaction between these clouds depends on the sign 
and size of Am = u x — u 2 , and on the product of the particle numbers MiM 2 - In our model, clouds are defined by 
the particles in a particular cell of a cubic grid. Two adjacent cells; which we will call a double-cell, will be picked at 
random and all the particles will undergo a collective collision which conserves energy and momentum in the double 
cell with a collision probability proportional to 8(Au) f (Au, M1M2) . It turns out that the equation of state can be 
tuned by the choice of the function /. 

Keeping the underlying lattice as a simple means to determine collision partners allows us to utilize the projec- 
tion operator formalism worked out for original SRD. This formalism was used to obtain transport coefficients via 
Green-Kubo relations Very accurate analytical expressions were found even for low particle density and at low 
temperatures 0, \2?\ where the mean-free path is much smaller than the cell size. Both regimes are usually avoided 
in DSMC, mainly because no good agreement between Enskog-theory and numerics can be obtained. 

II. MODEL 

Consider a set on N point-particles with continuous positions and velocities v^, i — 1...N. The particle mass is 
set equal to one. The discrete dynamics with time step t consists of a streaming and a collision step. In the streaming 
step, particles are advected freely with their corresponding velocities, ri{t + r) = Ti(t) + TVi(i). As in original SRD, 
0, E| a cubic grid with lattice constant a is imposed on the system, which is used to organize multi-particle collisions. 
In the current model, another grid with twice the mesh size of the original grid is introduced which groups four original 
cells into one supercell (see Fig. For simplicity, we restricted ourselves to two dimensions, but everything can be 
easily extended to more dimensions. Each supercell is indexed by the index of the small cell in its lower left corner. 
The origin of the grid is chosen such that the indices of the supercells are always odd. The position of a cell relative 
to its supercell is now uniquely determined by whether the components of its index vector ((, x ,£,y) are even or odd. 
For instance, the cell in the lower right has an even x-index, £ x , but an odd y-index, £ y . 

As proposed in 1] in order to avoid anomalies, all particles are shifted by the same random vector with components 
in the interval [—a, a] before the collision step (note this is a larger interval than in original SRD because of the bigger 
size of the supercells) . Particles are then shifted back by the same amount after the collision. 

To initiate a collision, two cells in every supercell are randomly selected. These two cells form what we will call 
a double cell. The probabilities for the possible double cells have to be determined in a way which is as efficient as 
possible, i.e. leads to a large non-ideal part in the equation of state, and furthermore should lead to a model which 
is isotropic at least at the Navier-Stokes level. As will be discussed later, the cubic anisotropy of the underlying grid 
will show up at the Burnett and higher levels, eventually leading to a cubic phase at very high collision frequency. 

For efficiency, we always pick two double cells in every supercell as possible candidates for momentum exchange. 
This way, every particle in the entire system has a chance to be included in a collision. As pictured in Fig. ^three 
distinct choices are made in every supercell: a) horizontal collisions: The two top cells with equal y-coordinate form 
one double cell, and the two lower cells form another double cell, b) vertical collisions: The two left cells with equal 
x-coordinate form one double cell, and the two right cells form another double cell, c) diagonal and off-diagonal 
collisions: The cell from the upper left and the one from the lower right form a double cell, while the remaining two 
cells form the second double cell. 

Note, that the diagonal collisions are essential to obtain an equilibration between the kinetic energies in the x and 
y-direction. Because of x-y symmetry, the probabilities for choice a) and b) must be equal, and will be denoted by 
w. The probability for choice c) is given by u>d = 1 — 2w. The ratio w/uid has to be determined such that the model 
becomes as isotropic as possible. One way to fix this ratio is based on the relaxation of the velocity distribution, and 
will be discussed in the next section. 
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FIG. 1: Collision rules. Three distinct collisions are considered: (a) horizontally along o\, (b) vertically along <T2, (c) diagonally 
and off-diagonally along cr-j, and CT4. w and Wd are the probabilities of choosing collisions (a), (b), and (c), respectively. 



A. Collision rules 

Fig. ^shows the three choices for double cells. We define unit vectors <Tj, j — 1..4, which connect the center of the 
selected cells. For the two possible double cells aligned with the x-direction one has cr\ = x and for the vertical choices 
one has <r 2 = y. The diagonal and off-diagonal choices are described by <x 3 = (x + y)/\/2 and <r 4 = (x — y)/y/2, 
respectively. In every cell the mean particle velocity is defined by 

1 m„ 

i— 1 

where the sum runs over all particles contained in a given cell, n is the cell index. 

Now, the projection of the difference of the mean velocities in the selected cells on <tj is calculated, Am = crj • (m — U2 ) . 
If Au < 0, no collision will be performed. For positive Au a collision will occur with an acceptance probability which 
depends on Au and the actual particle numbers in the two cells M\ and Mi- In order to enable analytical calculations 
and for first tests, the following total acceptance probability was chosen: 

PA = Q(Au) tanh(A) with A = AA11M1M2 (2) 

where A is some parameter to tune the equation of state. In the limit of A — > 00 this gives pa = O(Au), which 
maximizes the collision frequency, leads to a large non-ideal part of the pressure, and thus to a large speed of sound 
(the more collisions the faster sound can travel). Unfortunately, as will be shown in a later section, in this limit, 
the pressure has a non-analytic dependence on density and temperature, which leads to certain thermodynamic 
inconsistencies. More details about this limit can be found in Ref. 15]. Another, computationally simpler choice 
than Eq. would be the first term of an expansion of the hyperbolic tangent for small A: pa = min(l, 0(Au)A), 
but the cusp at pa = 1 makes accurate analytical calculations of the equation of state difficult. 
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This procedure represents a coarse-grained version of the collisions in a real gas. There, mostly binary collisions 
occur, and they are only possible if the particles approach each other, i.e. if Am = er 12 • (vi — v 2 ) > 0. v, is the velocity 
of the particles, <x 12 is the vector connecting their center of mass. In SRD the generalization of this can be seen as if 
two clouds of particles coming from the two cells are colliding. The total momentum transfer of this scattering should 
be much larger if the two clouds approach each other on average, i.e. at Au > 0, compared to Au < 0. Furthermore, 
the effective cross section of the scattering should increase with the particle density in the clouds. This is described 
by the dependency of the function A on the particle numbers. In this model there is large flexibility about how to 
choose this function with the restriction that it should be symmetric against the interchanging of the two cells. 

Once it is decided to perform a collision, the explicit form of the momentum transfer between the two cells is 
needed. In close analogy to the hard-sphere liquid, the collision should keep the total momentum and total kinetic 
energy of the double cell invariant, and it should mainly transfer the component of the momentum, which is parallel 
to the connecting vector crj. In the following, this component will be called the parallel or longitudinal momentum, as 
opposed to the perpendicular or transverse momentum. There are many different rules which fulfill these conditions. 
For instance, a stochastic rotation of the relative velocities of all particles in the double cells similiar to the rotation 
rules in SRD for ideal gases, could be used. Our goal here is to obtain a large speed of sound. Therefore, we use 
a collision rule which leads to the maximum transfer of the parallel momentum, and does not change the transverse 
momentum. The rule is quite simple, it exchanges the parallel component of the mean velocities of the two cells 
(parallel to crj), which is equivalent to a "reflection" of the relative velocities: 

^'(t + Tj-uN = -(«]'(*) -«") (3) 

it I is the parallel component of the mean velocity of the particles of both cells in a double cell. Written in Cartesian 

coordinates this amounts to: 

a) Horizontal double cells (characterized by cr-Cj: 



b) Vertical double cells (<r 2 ) 



c) Diagonal double cell (0-3): 



d) Off-diagonal double cell (0-4): 



v ix {t + t) = 2u x - v ix (t) 

V ly {t + T) = ViyQt) (4) 



v ix (t + T) = v lx (t) 

v iy {t + T) = 2u y -v iy (t) (5) 



Vix(t + r) = u x +u y -v iy (t) 

Vi y (t + r) = u x + u y -v ix (t) (6) 



Vi X {t + r) = u x -Uy + v iy {i) 

iHy(t + r) = u y -u x + v ix (t) (7) 

where u = (M1U1 + M 2 u 2 )/(Mi + M 2 ) is the mean velocity of the double cell, Ui, u 2 are the mean velocities of the 
cells forming the double cell. One easily verifies that these rules conserve momentum and energy in the double cell. 



B. Collision probabilities and isotropy 

The ratio w/wd describes how often collision cells in the vertical and horizontal direction, a) and b) are chosen 
compared to diagonal cell pairs, case c) in Fig. ^ We determined this ratio by requiring that the relaxation of the 
velocity distribution functions of the particles is isotropic and does not depend on the orientation of the underlying 
grid. 

Instead of analyzing the temporal evolution of the entire distribution function, we will restrict ourselves to its lowest 
moments. Assuming molecular chaos for the moment, i.e. that velocities of different particles are uncorrelated at 
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equal times, it is sufficient to consider the following three moments of a single particle i: (v 2 x ), (vf y ), and (vi X v iy ). For 
simplicity we will omit the particle label i from now on. We assume a homogeneous, but non-equilibrium inital state 
and study the temporal evolution of these moments to their corresponding values given by the Maxwcll-Boltzmann 
distribution. We define the vector of second moments as follows: 

/ («2(*)> \ 

^(t) = (V 2 y(t)) . (8) 

V M*H(f)) / 

For simplicity only the limit A — > oo in Eq. will be discussed in the following. This means, on average every 
second collision attempt is accepted (since the probability to have a positive velocity difference Au is 1/2). Using the 
collision rules, Eq. we obtain the following relaxation in one time step: 

*(t + r) =RV(t) (9) 

with the relaxation matrix R, 

/2w+(w d /2)(l + l/n) K/2)(l-l/n) \ 

R=\ (w d /2)(l-l/n) 2w+(w d /2)(l + l/n) . (10) 

\ w d + 2w/nj 

Here, n is the actual, fluctuating number of particles in a given double cell. Now we consider a rotation of the 
coordinate system by an arbitrary angle (3, i.e. 

* = ( ° 3 I * 

-s c 

with c = cos j3 and s = sin/?. This results in a transformation of the relaxation matrix from R to R = O R O^ 1 with 
a matrix 

2 2cs 

0=| s 2 c 2 -2cs ) . (11) 
-cs cs c 2 — s 2 

One of the three eigenvalues of R is always one due to the conservation of total energy, i.e. (v 2 ) + (v 2 ) = const = 
2ksT /m. The other two eigenvalues are given by 

2w 

Ai = w d H 

n 

A 2 = 2w+— (12) 
n 

Requiring isotropy amounts to R — R, which is possible only if Ai = A2, since the rotation of the coordinate system 
will mix these two modes. This condition can be fulfilled for arbitrary n only if w d — 1/2 and w = 1/4. Inserting 
these values we indeed find that the relaxation matrix is invariant under coordinate rotations: R — R. Another way 
to see this is that the requirement of R — R is equivalent to demand that R and O commute: [R, O] = RO — OR = 
This is only possible if R has the following structure with three free parameters ai, a 2 and 03: 

01 a 2 

R = I a 2 01 -a 3 I . (13) 



-a 3 /2 a 3 /2 a 




It turns out that original SRD with a rotation angle a is naturally described by a relaxation matrix of this shape and 
that the eigenvalues Ai and A2 are equal. In the current model, the probabilities w and w d must be adjusted properly 
to maintain this isotropic behavior. Adjusting these probabilities does not neccessarily mean that all properties of 
the model are isotropic. This becomes very apparent at high densities or high collision frequency 1/t>1 where one 
observes inhomogeneuous states with cubic order. Other practical tests for isotropy are needed such as measurements 
of the speed of sound and of transport coefficients for different directions of the wave vector k = (k x , k y ). Fig. [21 shows 
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FIG. 2: Adiabatic speed of sound c a measured using the time dependent density-density correlations for A/2a = 0.05. Densities 
are obtained at the supercell level and the period of oscillations in the density correlations are used to calculate the speed of 
sound. Data is shown for the wave vectors k = (2-7r/L)(l, 0) (•), k = (2tt/L)(0, 1) (□) and k = (27r/L)(l, 1) (a), respectively. 
Parameters: L/a = 128, A — > oo and ksT — 1.0. 



the speed of sound measured for three different directions of the wave vector k = (2-k / L)(n\, na), with (m, n^) — (0, 1), 
(1, 0), and (1, 1) where L is the linear dimension of the system. No dependence on the direction was detected, even 
at this small mean free path A/(2a) = 0.05, A = T^/ksT. A similar measurement (not shown here) revealed that the 
viscosity for a finite wave vector does not depend on the direction of k either. 

Choosing w = 1/4 and wa = 1/2 in Eq. 112(1 and averaging over the particle number fluctuations (particles are 
supposed to be Poisson-distributed as in an ideal gas) gives the effective decay rate, 



1 , ( 1 1 - e 



-2M 



Both, the cross-correlation (v x {t)v y (t)) and the difference between kinetic energy in x- and y-direction relax indepen- 
dently with this relaxation rate, i.e. 

(v x (t)v y (t)) = (v x (0)v y (0))e~ X e& f 
(v 2 x (t))-(v 2 y (t)) = ((^(0))-(^(0)))e- A eff t (15) 

Fig. [3] shows the time relaxation of these moments for small and large mean free path. The initial configuration was 
set up far from equilibrium with a large difference between (v£) and (Vy) and a strong correlation between the x- and 
y-component of the velocity. The measured decay is exponential for large mean free path A/(2a) = 1.5. For smaller 
A/ (2a) = 0.05 one sees a slower decay starting after a few iterations. The reason is that the assumption of molecular 
chaos is invalid at small A/a. This is similar to the behavior of the velocity auto-correlation function analyzed in |l6j . 

Fig0]shows the dimensionless decay rate ^ g qt as a function of M. At A/(2a) = 1.5 very good agreement is found 
for small and large M. At smaller A the decay is slightly slower due to the mentioned correlation effects. 
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FIG. 3: Non-equilibrium relaxation of velocity moments as a function of time, (a) Large mean free path A/(2a) = 1.50, (b) 
small mean free path A/(2a) = 0.05. Open circles (o) and squares (□) show the decay of (v x (t)v y (t)) and (u^) — (Vy) correlations 
respectively. The solid lines are the predictions of Eqs. ()14fl and 1151 . Parameters: L/a = 32, M = 5, A — » oo and fcsT = 1.0. 
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FIG. 4: Dimensionless relaxation rate A e g-r for the non-equilibrium relaxation of the velocity moments (v x v y ) and (v^.) — (vy) 
as a function of M, the mean particle number per cell, (a) X/2a = 1.50, (b) \/2a = 0.05. Bullets (•) show measurements and 
the solid line is a plot of Eq. I|14|l . Parameters: L/a = 32, M = 5, A — > oo and /c_bT = 1.0. 
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C. Phase-space contraction and detailed balance 

It is easy to see that the original SRD method conserves phase space and obeys detailed balance 0, 0|. The 
situation is more tricky for any coarse-grained, discrete-time method of hard-sphere collisions. In our model, any 
state in the 2dN-dimcnsional phase space is allowed, in contrast to a real hard-sphere system where particles cannot 
overlap. It is the collision rules described above which make certain states favorable than others. These rules ought 
to be biased to obtain a non-ideal equation of state and to build up non-ideal correlations between particle positions 
and their velocities. On the other hand, kinetic energy is conserved like in a real hard-sphere liquid, and we are in 
the microcanonical ensemble. Under these circumstances it is unavoidable that the dynamics must be phase space 
contracting, it has to project out certain states or at least sample them with a smaller probability. These unfavourable 
states correspond to the ones which are entirely forbidden in a real hard-sphere system. 

A superficial analysis seems to indicate that phase space is conserved: The Jacobian of the streaming step is always 
one and the Jacobian of the collision operation is either +1 or —1. A closer look reveals, that the dependence of the 
collision probability on Au destroyes this picture. Assume a N-particle state A with only one double cell and with 
Au < 0. No collision is allowed to happen, the state remains unchanged: A(t + r) = A(t) (We leave out streaming 
for the moment). Consider now a state B with identical particle positions as A but inversed velocities vf = —vf. 
This time Au > and a collision happens resulting in B(t + r) = A(t). Hence two different states are mapped onto 
exactly the same state in one iteration: phase space is contracted. 

As a consequence, the dynamics is not time-reversible even in a stochastic sense: Inverting the signs of all velocities 
of state B(t + r) = A(t + r) and performing one iteration would lead to a collision going back to the time-inverted 
version of state B even if one would have started at A initially. In other words, there is no way to go back in time to 
a state which did not fulfill the collision criteria, namely Am > 0. 

The whole dynamics presumably still fulfills detailed balance or at least semi-detailed balance, which is hard to prove 
since this involves not only the transition probabilities between states but also the (difficult to obtain) probability of 
the equilibrium states itself. Instead, we checked several thermodynamic properties of the system such as consistency 
of the thermal fluctuations with the pressure, [lil Ilfij. and the speed of sound which are related by thermodynamic 
expressions, and the scaling of interface fluctuations for the binary version of the model [iij. No inconsistency due to 
the absence of time-reversibility or due to a possible violation of detailed balance could be observed. 

Note, that also neither the consistent Boltzmann algorithm (CBA) |l7|. nor the coarse-grained model by Pooley 
j 2 ■ lt | nor any other method based on the Boltzmann or Enskog equation are time-reversible, for similiar reasons. 



D. Additional SRD step 

In order to change the viscosity independently from the equation of state, additional conventional stochastic rota- 
tions on the cell level can be performed. They will be denoted by the operator R. The streaming step with a time 
step r will be denoted by the operator S(t) and the biased collisions are called C. 

To obtain a time-reversal order of operations one has to be careful once the time evolution consists of more than 
two non-commuting operations. Consider the set of operations ... S C R S C R.... Inverting time, i.e. reading the 
sequence from the back leads to the chain ...S RC S RC... which is not equivalent to the previous one since C comes 
always after R, instead of R always following C. This set of operations can be symmetrized. One possibility would be 
using streaming with half the time step r/2 such as ...S(t/2) RS(t/2) C S(t/2).... This chain reads the same from 
the back, since R and C are embedded in S'-operations. Another computationally more convenient choice would be 
random symmetrization: ...SRC RS RC R.... Here R represents a rotation R which is only chosen with a probability 
1/2. This is the symmetrization we actually used in all our simulations. These ways of approximating time evolution 
operators are very similiar to Trotter- Suzuki formulaes used in modern symplectic Molecular Dynamics and Quantum 
Monte Carlo methods El. 
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III. EQUILIBRIUM CALCULATION OF TRANSPORT COEFFICIENTS 
A. Projection operator formalism 

In Ref. 0] a projection operator formalism was developed to derive the linearized hydrodynamic equations from 
the microscopic collision rules of Stochastic Rotation Dynamics (SRD). This technique was originally introduced by 
R. Zwanzig [23, [2II and later adapted for lattice-gases by Dufty and Ernst • With the help of this formalism, 
explicit expressions for both the reversible (Euler) as well as dissipative terms of the long-time, large-length-scale 
hydrodynamics equations for the coarse-grained hydrodynamic variables were derived. In addition, Green-Kubo 
(GK) relations were obtained which enable explicit calculations of the transport coefficients of the fluid. The GK- 
relations of SRD differ from the well-known continuous versions due to the discrete time-dynamics and the underlying 
lattice-structure. 

In the following, we will briefly outline how this technique is extended to accomodate the new collision rule, Eq. 
@, which in contrast to SRD does not conserve momentum and energy in single cells but in randomly chosen double 
cells instead. For more details about the formalism the reader is referred to Ref. 

The starting point of this theory are microscopic definitions of local, hydrodynamic variables Ap. These variables 
are the local density, momentum, and energy density. They can be defined on the cell level as 

N d 

Mi) = E a ^ IT (I - & " r*rl) > ( 16 ) 

i=l 7 = 1 

with the discrete cell coordinates £ = am, with mp = 1,...,L, for each spatial component, an — 1 describes 
the particle density, {ap.i} = {"Vi(p-i)}, with (3 — 2,...,d + 1, are the components of the particle momenta, and 
a d+2,% — v i/% is the kinetic energy of particle i. d is the spatial dimension, r, and Vj are position and velocity of 
particle i, respectively. The hydrodynamic variables look much simpler and are easier to treat in Fourier-space. Their 
spatial Fourier transforms are given by 

^(k)=^a w eH, (17) 

3 

where £ ■ is the coordinate of the cell occupied by particle j. k = 27rn/(aL) is the wave vector, where np — 
0, ±1, . . . , ±(L — 1), L for the spatial components. 

The next step is to set up an evolution equation in discrete time which looks like a continuity equation in Fourier- 
space, 

A t Ap(t)+ik-T)p(t) = , (18) 

where A t is the discrete time evolution operator, defined as A t A(t) = [A(t+T)—A(t)]/T. However, up to this point this 
equation does not contain any information; the flux Dp is just formally defined by ik- Dp(t) = — [Ap{t + r) — Ap(t)]/T. 
The key point in making this expression meaningful is to set up a microscopic expression for the local conservation 
of density, momentum and energy (for SRD this expression is given by Eq. (23) in Ref. 0). This expression is then 
inserted into the formal definition of the fluxes in order to cancel divergent terms proportional to 1/k. If this can 
be done, the reformulated Dp really is a flux and Eq. i|18|) is proven to be a continuity equation for the conserved 
quantities Ap. 

In the current model, for every choice of double cells, the local conservation law can be expressed as 

(e ik -«K*+^) + e ik-(«J(t+r)+»,«)) [o w (t + r)-o w (t)] = 0, (19) 

3 

where is the coordinate of the cell occupied by particle j in the shifted system. Note, that similiar to original SRD, 
a random shift of cells before collisions is required to remove anomalies at small mean free path, see Ref. 0. The 
vector Zji is a function of the cell coordinate £j and has components which are either 0, 1 or —1. It is constructed 
such that the sum of the two exponentials in Eq. I|19|) is the same for two particles if and only if they are in the 
same double cell. The index 2, I = 1, 2. ..6 describes the choice of double cells in more detail than the vector rrf I = 1 
denotes the lower horizontal double cell, whereas I = 2 describes the top horizontal double cell in Fig. ^ Similarly, 
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/ = 3 stands for the left vertical double cell, and I = 4 is needed if the particle j happens to be in the right vertical 
double cell. I = 5 is for the diagonal choice, and I = 6 for the off-diagonal choice of double cell. 
The components of Zji are defined by: 




odd, 
even, 



even 
odd 
odd 



(20) 



Zj2,x — 



odd, 
even, 



&v odd 

e 



■s 



even 
even 



(21) 
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even, 
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odd, 




odd 
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odd, 
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even 
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( I 




odd, 






z 3^,y 
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even, 




odd 








£ s 


even, 




even 
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£ s 


even, 




odd 











odd, 


^3,x 


even 


^jb,x 




1 




odd, 




odd 






\-i 




even, 




even 




= 2 
















( 




even, 




even 


Zj6,x 









odd, 




odd 




1 




odd, 




even 










even, 




odd 



(22) 



(23) 



(24) 



(25) 



In summary, Eq. I|19fl is just another way of saying that the mean particle number, the mean velocity and energy in 
a double cell is the same before and after a collision. 

For the moment we will assume that only one choice for the construction of double cells within a supercell is possible, 
i.e. I is kept fixed. Later, an average over I with the probabilities w and Wd will be performed. To obtain the evolution 
equations for the hydrodynamic variables Ap in k-space, the differential quotient A t Ap = [Ap(t + t) — Ap(t)]/r has 
to be written in a form ~ ikD^, where the flux Dp is of order 0(k , ). Hence we have 



rk D, 



Ap(t + r)-Ap(t) = 2^ \ap d e 



(26) 



with apj — ap,j{t + t). The first term on the r.h.s. can be written as: 



ik-e 3 (t+T) 



3 



3 ik-(£f (t+r)+z ji ) + e ik-(« 3 s (t+r)) 



el k-(s 3 s (t+r)+z 3i ) _ e tk-(«f (t+r)) 



(27) 
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where we have added and subtracted two identical terms. A similiar identity can be applied to the second term on 
the r.h.s. of Eq. 



j 



a /3li (t)e <k -«'W = ^ap tj (t) (V k ^ (t) + - \ e ^Uf(t+r)+^) + „&.(«?(*+,-)) 
j 2 



ik-( Sj s ( t+r ) +z . i ) _ ik-(«J(t+r)) 



(28) 



As before, the quantity in the square brackets is identical to zero. Now, both expressions H27|) and i|28|) are inserted 
back into Eq. (|26|1 . Terms coming from the second and third terms in the r.h.s. of Eqs. I|27|) and (|28|) cancel each 
other because of the local conservation law, Eq. I|19|l . The remaining terms can be rearranged in the following way, 

1 
2 



A fi (t + r)-A p (t) = E(««- a ^) elk ' Sj(t+T) 



e ik-(tf(t+T)+x jt ) + e ik.(«f (t+i-)) 



a ^(*) ( ( 



(29) 



This expression is of order k for small k. Comparison with the definition of the flux Dp in Eq. (|26[1 and expanding the 
exponentials around £j(t) for small k shows that the flux Dp is indeed of order 0(k°) as expected for hydrodynamic 
modes, 



D/> = ~E {o W (t)^(t) + Aa w (t) [A«5(t) - ^] + O(zk)} e^W 



(30) 



where A£ ■ = ^ •(< + r) — A£ ■ = ^ -(i + r) — £ ■ (t + r), and Aap t j = ap,j{t + t) — aa,j(t). This expression turns 



out to be a simple extension of the flux of the original SRD fluid (Eq. (50) in part I the only difference being 
the vector Zji which describes the coupling to a neighbor cell. Hence, the entire formalism for calculating linearized 
hydrodynamics and transport coefficients for an ideal SRD-fluid can be applied without changes except the modified 
definition of the flux. In particular, using Eq. (66) of [6j one easily obtains a Green-Kubo relation for the kinematic 
shear viscosity v. 

Using the resummation procedure described in 0, this GK-relation can be rewritten in a form which is easier to 
evaluate analytically. The resummation utilizes the exact Galilean-invariance provided by using random grid shifts 
and cancels all terms containing the unshifted cell positions £j. Then, the viscosity is obtained in terms of the 
off-diagonal part of a resummed stress tensor a xy , 



Nk * T h 



(0)a X y(nT)) 



The prime means that the first term in the sum is multiplied by a weight 1/2, and 



N 

Oxyinr) = ^2 \ v jx(n>i~)v jy (nT) + v 3y (nT)B jx (nT)} 



(31) 



(32) 



3=1 



where 



B 3X {nr) = ef x ((n + l)r) - ^(nr) - r^(nr) + ~(z jx ((n + l)r) - z jx (nr)) 



(33) 



Note, that the only formal difference to original SRD jj] is the appearance of the z variables. Therefore, we write 



B(nr) = B(nr) + -[z((n + l)r) - z(nr)] 



(34) 



where B(nr) is the definition used earlier in Ref. |4|. 

An important test for the validity of the new model and the definition of the flux is, whether the reduced frequency 
matrix u> (see Eq. (48) in 6]) is symmetric. In we checked that uo is symmetric for the original model. This 
matrix contains scalar products of the form (Ap\Dry), (D^A^) in the limit of k — > 0, i.e. averages where the flux and 
hence the vector Zji occur only linearly. Note, that Zji is uncorrelated to the velocity- variables ap t j- Furthermore, the 
probability is 1/2 that £ ja is even or odd. Hence, (zj/) = for all I, and (zjiap, n ) — 0. It turns out that the additional 
non-ideal term in the flux does not contribute to the reduced frequency matrix; the matrix is symmetric and identical 
to the one of the ideal SRD-fluid. 
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FIG. 5: Viscosity as a function of time step r. Filled circles (•), squares (■) and triangles (A) show results for M = 1, 3 and 5 
respectively. Measurements are done fitting vorticity decay profiles for the smallest wave vectors and averaging over 5 different 
ensembles. The error bars reflect this averaging. The inset shows a comparison of these vorticity measurements (■) with the 
Green-Kubo measurements (□) for M = 3. The solid lines are plots of the sum Vkin + v co n given by Eqs. 13511 and 1371 in both 
graphs. Parameters: L/a — 32, A — > oo and ksT = 1.0. 



B. Determining the viscosity 

The Green-Kubo relation given by Eq. Il-'S1[1 was evaluated analytically and numerically. If only the first term in 
the stress tensor, Eq. (13211 . (which contains terms ~ Vj X Vj V ) is used in the evaluation of 1)31(1 one obtains the so-called 
kinetic viscosity Vkin ■ Including particle number fluctuations based on the Poisson-distribution, we found in the limit 
of large acceptance probability A — > oo (see Eq. (J2J): 

fc^Tr / 6M + l-o-™ \ 
v Un~ 2 I^M-l + e- 2 ^ [ ' 

M is the average particle number in a cell of size a. For small acceptance probability i.e. for small A, one finds 



^AM 3 / 2 \j k B T 2 

where we have neglected particle number fluctuations. Since the first term in Eq. I|36l) dominates at typical acceptance 
rates, the viscosity is proportional to the square root of temperature as in a real gas. Interestingly, in this limit it 
turns out that the self-diffusion constant is exactly equal to the kinetic viscosity, D = Vkin- More details about the 
derivation and numerical comparisons will be published elsewhere. 

For large mean free path X/ {2a), the total viscosity is just given by the kinetic viscosity. However, at small mean 
free path, additional contributions, denoted by v co ii become relevant. For A — > oo this contribution takes the form 

a 2 ( 1 - cr 2M \ 

*— M 1 — sH (37) 
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However, the prefactor G is not known yet and might depend on time step r and temperature T. Exploring these 
details is subject of ongoing work. 

Figure [S] shows measurements of the total viscosity as function of mean free path and mean particle number per 
cell. For comparison, the theoretical expressions for i/£° n + v co u from Eq. 1351) and 13 7|) are shown with an arbitrary 
constant prefactor G — 1/6. Good agreement is found for mean free paths A/(2a) > 0.25. The inset in Fig. 
compares measurements of the viscosity using two very different concepts, (i) using vorticity correlations as explained 
in Ref. 0, and (ii) using the GK-relation l|31l) . The results of both methods agree within a few percent. 



In two dimensions, the velocity-autocorrelation function R(t) = (vi X (0)vi X (t)), and the stress correlation function 
(cr X y(0)<Jxy(t)) are expected to decay as \/t at large times. Since the transport coefficients are related to the integrals 
of these correlation functions by GK-relations, one should see a logarithmic divergence of the time-dependent viscosity. 
These "long-time tails" and the logarithmic divergence have been measured in the original SRD-method and found 
to be in quantitative agreement with general hydrodynamic theory [in], [H, HH . However, as can be seen in Fig. 
no such divergence could be detected in the present model. Instead, the viscosity quickly converges to a plateau and 
does not drift with time. The amplitudes of the expected divergent terms were calculated analytically in the 70's 
[siL Is^ and shown to be essentially inversely proportional to the density and to the transport coefficients such as the 
shear viscosity and the self-diffusion constant. We estimated these amplitudes and found that they are too small to 
be seen within our numerical resolution, since the present model is more viscous than original SRD. 

Fig. Hshows the velocity auto-correlation function. For small A, a calculation based on the assumption of molecular 
chaos gives the following exponential behavior: 



The measured decay of the auto-correlation function R(t) is exponential over four decades and agrees perfectly with 
expression |J3S}, i.e. there is no sign of long-time tails. 



In order to obtain explicit hydrodynamic equations, the susceptibility matrix must be known. These equal time 
correlations do not follow from the projection operator technique and must be determined by other means. For a hard- 
sphere fluid, the only quantity needed to obtain all necessary information is the equation of state. The other quantity, 
the internal energy, is already known since it is identical to that of an ideal gas. This is because we constructed our 
collision operations to leave the kinetic energy invariant in a double cell. Here, we follow the mechanical route to 
pressure, which is defined as the average longitudinal momentum transfer across a fixed interface per unit time and 
unit surface. 

Let us assume this interface to be parallel to the y-direction, and only consider transfer of x-momcntum, i.e. the 
component p xx of the pressure tensor is to be determined. Double cells aligned with the y-direction do not contribute, 
only double cells characterized by <Ti, er 3 , and er^ have to be considered. We will only discuss momentum transfer 
caused by collisions, but not the contributions from streaming which give rise to the ideal part of the pressure. 

Since the collision rules are applied in homogeneously shifted cells, the distance S x between the left most corner of a 
double cell and the dividing line in x-direction is homogeneously distributed between and 2a. The average amount 
of momentum transferred across the dividing line is zero for 5 X = and S x — 2a, and increases linearly towards 5 X = a, 
i.e. reaches a maximum if the dividing line goes through the center of the double cell. Averaging over this position 
dependence leads to Gtot = ^G m id- Gtot is the averaged transferred momentum, G m id is the transferred momentum 
across a line going through the center of mass of the double cell. 



C. Absence of long-time tails 




(38) 



IV. THE EQUATION OF STATE 



A. The mechanical route to pressure 



10 





FIG. 6: Temporal evolution of the shear viscosity, (a) A/2a = 1.50, (b) A/2a = 0.05. The solid and dashed lines show the 
kinetic and total contributions to shear viscosity as a function of time step, respectively. The dotted lines show the prediction 
of Eq. for "iSn- Parameters: L/a = 32, M = 5, A -> oo and k B T = 1.0. 
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100 



FIG. 7: Normalized velocity auto-correlation function as a function of time for A/2a = 0.25. The solid line is the prediction of 
Eq. Parameters: L/a = 32, A = 1/60, M = 5 and k B T = 1.0. 



Let us consider horizontal double cells first. The change in x- velocity of a particle in one cell during an accepted 
collision is given by Q), 



Avi X =2{u x - v ix ) 
hence, the x-momentum change in the left cell is 

AG X 



2M 1 M 2 
m l\u, 



M 1 +M 2 



(39) 



(40) 



Au = u\ x — u 2xi ui x and u 2x being the average x- velocity in the left and right cell, respectively, m is the mass of the 
particles which is set to one. 

For the moment we will assume, that the particle numbers Mi and M 2 will be constant, i.e. the following averages 
are performed under this restriction. To obtain the pressure, the thermal average over the momentum transfer across 
a line is required which involves knowledge of the acceptance probability pa{M\, M 2l Au) for a collision. One obtains 



W / 

(AG,) = - / PG {Au) p A (M 1 ,M 2 , Au)AG x d(Au) 
* Jo 



(41) 



The factor 1/2 comes from the position average of the dividing line; the integral involves only positive Au, because 
the acceptance rate is zero for Aw < 0. pg(Au) is the probability that u\ x — u 2x for the micro-state of two cells is 
equal to Au. w is the probability to select this double cell and will be equal to 1/4 as discussed before to ensure 
isotropy. One has 



Pg(Au) 



oo M1+M2 I ^ Mi M1+M2 

"°° »=1 \ 1 j=l Z 3=Mi+l 



(42) 
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M 

FIG. 8: Non-ideal part of pressure times r as a function of M. Data for r = 0.05,0.10,0.20,0.40,0.60,0.80,1.00,2.00 are 
collapsed. The solid line is the first term (~ fcsT) of the theoretical prediction of Eq. 11491 . Parameters: L/a = 32, A = 1/60 
and k B T = 1.0. 



where g(yi X ) is the Boltzmann- weight, g{vi X ) ~ e~xp(mv~ x /2ksT) . Using the integral representation of the 5-function, 



M.r) I e ifcc ^ (43) 



Eq. igU leads to 



p G (Au) = J—j——e^BT , (44) 



with 7 = (Mi + M 2 )/MiM 2 . 

Expanding the acceptance probability, Eq. in A = A Au M1M2 leads to 

PA (M l ,M 2 , Au) = 8(Au) U~y + -) ■ ( 45 ) 

The contributions to the pressure from all terms of this series can be exactly calculated but we restrict ourselves to the 
first two terms. The reason is that all terms which are non-linear in A are related to a thermodynamic inconsistency 
and have to be kept small anyways. The size of the first non-linear contribution will give us clues about the useful 
parameter range and estimations of the violation of thermodynamics. 

Defining the pressure as the average momentum transfer per unit area and unit time, p\ — (AG X ) j (ra d ^ 1 ) [d is 
the spatial dimension, which is two in our case), the pressure follows from Eqs. I|41I44|I as 

u)A/c B T r „ r wA 3 (k B T) 2 „,o„,o 
pi = — —Mi M 2 - -— y B ' 1 M\ Ml + ... (46) 
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The index 1 for the pressure denotes the contribution from horizontal double cells. A calculation following the same 
lines yields that p% (due to the diagonal collision <x 3 ) can be obtained by just replacing w by w,j/V2 in the expression 
for pi and by symmetry we have ps = p±. 

Using w = 1/4 and Wd — 1/2, the total non-ideal pressure under the restriction of the particle numbers in the two 
cells to be Mi and M 2 is therefore 



P 



( 1 



V2V2 



A k„T 



2 w 



d-l. 



Mi M 2 



A 3 (k B T) 



Ml) 



(47) 



For small particle densities M n < 4, the fluctuations of the particle number in a cell cannot be neclected. 
As a first approximation, we assume the particles to be Poisson-distributed, 



p{M 1 ) 



-M 



M Ml 
"MjT' 



Mi = 0, 1, 2. 



(48) 



p{Mi) is the probability to find Mi particles in a given cell, where M = (Mi) is the average number of particles in 
that cell, and we further assume the distributions of adjacent cells to be independent: p(M\,M2) = p{Mi)p(M2). 
These assumptions are only strictly true for an ideal gas, but turn out to be sufficiently accurate. In principle, the 
corrections of this to the non-ideal equation of state and the resulting correlations of the particle numbers in adjacent 
cells should be calculated self-consistently. 

With these approximations, one obtains the non-ideal part of the pressure averaged over particle number fluctuations 
in two dimensions: 



id 







A M 2 k B T 






2 ar 



1 - 2A 2 M 3 k R T 



4 
M 2 " 



1 

M3 



(49) 



Note, that P non -id is quadratic in the particle density p — M/a to lowest order in p and for small A as one would 
expect from a virial expansion. Furthermore, the first term in the non-ideal pressure is linear in temperature T which 
is also expected for a hard-sphere gas. The second term is just the first contribution arising from the non-linearity of 
the hyperbolic tangent in the acceptance probability. All these terms lead to a non-linear dependency on temperature, 
which will be shown to be inconsistent with the internal energy of a hard-sphere gas. This means, the prefactor A 
has to be chosen small enough to be in the linear regime. This is equivalent to having a small acceptance rate for a 
collision. It turns out that acceptance rates of about 15% are sufficiently small to avoid inconsistencies. 

In the limit A — > 00, the acceptance rate is 50% because Au has to be positive to accept a collision. In this limit and 
for large M an expression can be derived which shows a non-analytic dependence of both density and temperature, 
namely P n on-id ~ \Zk B Tp which again is thermodynamically inconsistent as discussed above. Another technical 
difficulty for this case is that the average over the particle fluctuations cannot be performed exactly anymore, only at 
large densities an asymptotic result is found. On the other hand, even for arbitrary values of A the pressure can be 
obtained as an infinite sum by expanding the hyperbolic tangent. All terms in this sum can be calculated exactly. 

The off-diagonal terms of the pressure tensor can be calculated along the same lines. One finds p\ xy — 0, pzxy = 
Pixx — —Pixy Hence the total off-diagonal term P xy vanishes. In the limit of A — > 00 and for large M, the equation 
of state has the following non-analytic form: 



poo 
non—id 



= pk B T 1 + - 



1 1 



A JM 



1 



1 



8 4V2 



(50) 



where a is the lattice unit of the cell, A is the mean free path A = T\/k B T. 



B. A microscopic expression for the pressure 

The microscopic expression for the flux D/3 was derived in Eq. (|30J) . The components /3 = 2, ..1 + d also provide a 
microscopic expression of the stress tensor: the average of the diagonal part gives the virial expression for the pressure, 



Ptot — Pideal 
1 

7v 



Pnon — i 



^ |» V., AC, + Ar, 



AS 



%X,jl 



}) 



(51) 
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The first term was discussed in part II of Ref. p, (v x ,jA^j) = (rv%. -) and gives the ideal part of the pressure, 
Pideai- The average over the second term vanishes (see Ref. and it is the third term leading to an expression for 
the non- ideal part of the pressure. As seen directly from Eq. I|51l) . this term counts the momentum transfer between 
the two cells in a double cell: 

Consider a horizontal double cell for simplicity. Let — Ap be the total momentum change (in x-direction) in the 
left cell. Due to momentum conservation in the double cell, the total momentum change in the right cell is equal to 
Ap. Now, z x is equal to 1 in the left cell and equal to —1 in the right cell. Hence, as a result of the multiplication by 
z x /2 in Eq. (|51l) we have Ap(l/2) — Ap(— 1/2) = Ap, i.e. the non-ideal part of the pressure is expressed as the sum 
of all the momenta which are transferred from the left cell to the right cell in all selected double cells in the system. 
This transferred momentum is mostly positive, since only the collisions with a higher mean velocity in the left cell are 
allowed to happen. This explains why biased collisions are needed to obtain a non-ideal equation of state; in original 
SRD there is still momentum transfer due to rotations across a plane but with positive and negative contributions 
which cancel each other. 

In the current model, due to particle number fluctuations, there can be situations where the mean velocity in the 
left cell is higher than in the right allowing a collision, but the mean momentum = u\ M\ is smaller than in the right 
cell with Mi < M 2 . This would lead to a negative contribution to the pressure. However, these are events of small 
probability and are even less relevant at higher particle density. On average, unlike in original SRD, the amount of 
transferred momentum does not vanish, and we have indeed a non-zero pressure. 

Comparison with the analytical expression, (|49|l shows good agreement, Fig. |S| There are deviations at large M 
because the acceptance rate gets larger with M, and at rates above 15% the terms proportional to (fc^T) 2 become 
important in the pressure, 149|) . These corrections were not included in the theory plotted in the figure because the 
model becomes thermodynamically inconsistent in that range. Better agreement at larger M can always be obtained 
by decreasing the prefactor A in Eq. J2J) . 



C. Density fluctuations and thermodynamic consistency 

Thermodynamics gives a relation between derivatives of the pressure p and the structure factor S(k, t) at zero time 
t = and in the small wave number |fc| — > limit, 



S(k,t = 0)= P k B T„ 
op 



(52) 

T 



On the other hand S(k,t = 0) is given by the equal time correlations of the density fluctuations in Fourier space, 
(pfe(O)p-fc(O)). Hence, a numerical consistency check can be performed: The density fluctuations are measured, and 
the derivative of the measured pressure is taken numerically and inserted into Eq. I|52l) . Both routes should lead to 
the same S(k, t = 0) which can also be calculated analytically using the equation of state, l|49|l . As shown in Fig. 
this works nicely for small prefactors A and both large and small time steps r. For small M (< 5) good agreement is 
achieved. For larger M, one comes closer to the limit of A — ► oo, where the acceptance probability is independent of 
particle densities; this limit is thermodynamically inconsistent as we will discuss below. Of course, if agreement for 
M > 4 is required, all one has to do is to reduce A, i.e. decrease the acceptance probability. 

As mentioned before, at large acceptance probability, where the non-linear terms in A in the equation of state are 
not negligible, the model is not thermodynamically consistent. The reason is that these terms contain non-linear 
functions of temperature ksT. This leads to the following problem: let us assume a term in the equation of state of 
the form T 2 p n , where 1. The pressure is related to the free energy density / by the relation 

P(T,p)=p& T -f . (53) 
dp 

This is a linear differential equation for / with the general solution 

f = C(T)p+—P- (54) 
n — 1 

where C is an arbitrary function of temperature. The entropy density is denned by the following derivative of the 
free energy density: 

s = ~%\p (55) 
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(a) _ 




FIG. 9: Static structure factor S(k,t = 0) as a function of (a) t for M = 3, (b) M for r = 1.0. The open (o) circles show 
measurements from taking the numerical derivative of the pressure. The filled (•) circles show direct measurements of the 
density fluctuations. The solid line shows the theoretical prediction based on Eqs. 1491 and 1521 . k denotes the lowest possible 
wave vector, k = (2n/L) (1, 0). Parameters: L/a = 32, A = 1/60 and fcsT = 1.0. 
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which gives in our case: 



dC 2T P n r«n 
ol n — 1 



Finally, we obtain the part of the internal energy density related to the pressure term ~ T as 

The first term is linear in p, the second one is nonlinear in p. Hence it is impossible to have u — at all densities even 
for an arbitrary function C(T). This is a contradiction, because we have a model with collision rules which exactly 
conserve kinetic energy. Hence, the internal energy should be the same as for an ideal gas, i.e., there should be a 
kinetic term only and the additional contribution related to the non-ideal part of the pressure should be zero at any 
temperature and density. The fact that there are additional contributions to the internal energy probably means that 
the thermodynamic temperature does not agree with the kinetic one (defined via the square of the particle velocities) 
and might be the cause of the observed inconsistencies in the density fluctuations at large A (defined in Eq. J2J). 
This seems to be supported by another contradiction we observed at large A: The measured temperature fluctuations 
seem to be consistent with the specific heat cy = dfcs/2 of an ideal gas, but the theoretical prediction for cy based 
on a free energy density with non- linear temperature dependence differs from dfcs/2. 
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FIG. 10: Ordering at small r. Initial configuration is 23 x 23 clouds. Observed configuration is 23 x 23. Parameters: L/a 
A = 1/60, M ~ 4.69, r = 0.0005 and k B T = 1.0. 
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V. CAGING AND ORDER/DISORDER TRANSITION 



When the non-ideal part of the pressure is significantly larger than the ideal part, ordering effects can be expected. 
As in a real hard-sphere gas, both parts scale with temperature in the same way (for small A). Hence, as in a real 
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system which does not have an energy scale, we can assume that changing the temperature alone will not lead to an 
order/disorder transition — unless of course temperature becomes so large that the nonlinear terms in the pressure 
become relevant. On the other hand, the two pressure terms have different dependencies on time step r and density. 
t can be interpreted as a parameter describing the efficiency of a collision; lowering r would have a similar effect as 
making the spheres larger in a real system, resulting in a higher collision frequency. Therefore, we expect caging and 
ordering effects if either p is increased or r is decreased. This is indeed the case. For r < 0.0016 and M w 5, an 
ordered cubic state is observed, as shown in Fig. ^| In Fig. ^2 the pair-correlation for this state is plotted. Unlike 
in a real system, clouds of several particles are concentrated at locations r nrn = na x + ma y , where the 
periodicity in x and y-direction, respectively. The average number of particles in these clouds will be called the cloud 
number. We found that there is usually at least four particles per cloud. These ordered structures are similar to the 
low-temperature phase of particles with a strong repulsion at intermediate distances, but a soft repulsion at short 
distances. 

One of the surprising features of this crystal-like state is, that x-y symmetry can be broken, i.e. a x is not always 
equal to a y . Furthermore, there is the possibility of having several possible metastable crystalline states corresponding 
to slightly different lattice constants and cloud number. For instance, Fig. 1 101 and Fig. ^1 were created for the same 
parameters but different initial conditions. While there is 23 x 23 clouds in the first figure, there is 20 x 19 — 20 in 
the other. This means also that in the latter picture a x and a y differ slightly. Both states are relatively stable over 
long times. However, in Fig. ^|we see a lattice defect (circled for better visibility) and there is still activity in the 
lower middle part of the sample where particles seem to migrate from one cloud to another. As expected, the lattice 
constants a x ,a y « 1.6 are slightly smaller than the super cell spacing 2a which sets the range of the multi-particle 
interaction. In this state the diffusion coefficient becomes very small, particles are caged and can barely leave their 
location. As shown in Fig. 1141 we measured the coefficient of self-diffusion D as a function of time step r in the 
long-time limit, i.e. after an eventual transition to an ordered state. Above a critical tc ~ 0.0016 (at M = 5) the 
measured value of D agrees perfectly with the theoretical prediction: 
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FIG. 11: Pair-correlation function for the configuration shown in Fig. 1101 
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FIG. 12: Ordered state with defect. Initial ordered configuration is 16 x 16 clouds. Observed configuration is 20 x 19 — 20. 
Parameters: L/a = 32, A = 1/60, M ~ 4.69, r = 0.0005 and k B T = 1.0. 




FIG. 13: Pair-correlation function for the configuration shown in Fig. 1121 
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FIG. 14: Self-diffusion coefficient D as a function of time step r. The solid line is the prediction of Eq. I|58|l . Parameters: 
L/a = 64, A = 1/60, M = 5 and fc s T = 1.0. 



(A is small and particle number fluctuations were neglected here). However, D dramatically drops below tq which 
is exactly the value of r below which we clearly see an ordered state. This pronounced effect allows us to locate 
the transition point quantitatively. Another way to describe the transition is by determining the pair-correlation 
function g(r). Fig. 1151 shows that the behavior of g(r) is very different slightly above and below tc- Below tc we see 
oscillations due to a cubic structure with long-range order with lattice constant « 1.6 a. To describe the transition in 
more detail one could plot the amplitude of the Fourier-mode for the oscillation in g(r) (not shown here). 

How can we understand the stability of the ordered state? Without collisional interaction, particle clouds will 
broaden due to streaming, which will happen the faster the higher temperature is. Due to the grid shift, particles 
at the perimeter of the clouds will more often undergo collisions with neighboring clouds. These collisions provide 
back-scattering, particles are forced to fly back towards the center of their cloud. Thus, a correlation between distance 
from cloud center and speed is built up leading to stable cloud formation. 

VI. GENERALIZATION TO A BINARY MIXTURE 

The collision model can be easily generalized to multi-component mixtures. Consider a binary system with two 
types of particles, A and B. In order to obtain phase separation we introduce a repulsive interaction between different 
kind of the particles, but no repulsion among particles of the same kind. This is done in the following way: Suppose 
a double cell is selected for a possible collision. A-particles in cell 1 can now undergo a collision with B particles from 
cell 2. Furthermore, B particles from cell 1 are checked for possible collision with A particles in cell 1. The rules and 
probabilities for these collisions are exactly the same as in the one-component situation. Since in a phase separated 
situation there is hardly any collision away from the phase boundary, additional regular SRD-rotations on the cell 
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level are incorporated to mix particle momenta. 

A phase separation into an A and B rich phase is observed above a critical prefactor A of about A c = 0.36 for 
Ma = Mb = 5, Ma.b being the average number of A and B particles in a cell, respectively. The phase diagram does 
not depend on temperature as expected since similiar to a hard-sphere gas there is no energy scale in the present 
model. The spectrum of the interface fluctuations of droplets well above A c was analyzed, and it was found that 
they scale with wave number k as 1/fc 2 as expected. Quantities like pressure, chemical potentials, surface tension and 
transport coefficients were also obtained analytically and will be published elsewhere [25| . 



VII. CONCLUSION 



In this paper we analyzed a particle-based model for fluid dynamics with effective excluded volume interactions 
which was introduced in Ref. |b|. These interactions are modeled by means of stochastic multiparticle collisions 
which are biased and depend on local velocities and densities, but exactly conserve momentum and energy locally. The 
motivation for such a model and the collision rules were explained in detail. The relaxation to thermal equilibrium 
was investigated analytically and numerically, and it was shown how this relaxation was used as a guide to make the 
model as isotropic as possible on a cubic grid. A brief summary of how to calculate the pressure was already given 
in Ref. 0|; here we give a more detailed description of the calculation. We also show how a microscopic formula for 
the pressure can be derived, which allows local measurements of the pressure. 

In addition, we outlined how a discrete-time projection operator technique developed for original SRD 6j can 
be generalized to the current model, and how Green-Kubo relations can be obtained for the transport coefficients. 
Numerical measurements of the velocity and stress auto-correlation functions did not show long-time tails even though 
all the simulations were done in two dimensions. We think that this is mainly due to the large viscosity of the model, 
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leading to amplitudes of the tails much smaller than the numerical resolution. 

In the limit of very large collision frequency, i.e. at large particle density and/or small time step r we found an 
ordered cubic phase. The four-fold symmetry of this state is a numerical artifact due to the underlying cubic grid, 
but the crystallization itself resembles caging in a real gas and is a generic feature of models with soft repulsion at 
short distances. This order/disorder transition was quantitatively located by measuring the self-diffusion coefficient 
D. At the transition we saw a pronounced change of the behavior of the diffusion coefficient: D is about one to two 
orders of magnitude lower in the ordered phase than it would be in a homogeneous system. 

The ordered state consists of particle clouds containing at least four particles. These clouds form metastable cubic 
phases with long-range order. The lattice constant of these phases could slightly vary in x- and y-direction for the 
same set of parameters. This is due to the additional degree of freedom of how many particles on average can live in a 
cloud. Lattice defects were observed leading to a slow evolution from less favourable states where the lattice constant 
was below or above the optimum range. This range seems to be centered around 1.6a. 

The current model was extended to multi-component mixtures such as binary mixtures and microemulsions, and 
the critical behavior was investigated. These studies will be published elsewhere. 
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APPENDIX A: AVERAGES OF Au 



Averages of Au, the mean velocity difference between two cells, show up in the calculation of the acceptance rate 
for a collision. They can be measured as a test of the correct implementation of the code. We restrict ourselves to 
horizontal cells, where only the x-componcnt of velocities are needed. In this case Au is defined as 

a Vix + V 2x + ••• + V Ml + l,x + V Ml +2,x + ■■ + V Ml + M 2 ,x /A1 s 

Au = — — '■ — . (All 

Mi M 2 

Here M\ and M 2 are the current particle numbers in the two involved cells, respectively. The average over the velocity 
distribution at fixed M\ and M 2 is given by 



(|Au|) = (|°°A« p G (Au) dAu) = y ^ (A2) 

where pa is given in Eq. I|42(l . The particle mass is set to one as before, and 7 is defined as in Eq. Ij44(l . 

The average over the particle fluctuations, which we again assume to be Poisson-distributed, and uncorrelated on 
the cell level cannot be done analytically. Let us first consider ((Am) 2 ) which does not have the complication of a 
square root, y/j. It involves the following averages 

n— 1 

The double bracket denotes the average over the particle number in a cell, while the single bracket is for averages over 
velocity distributions. The r.h.s. of Eq. l|A3Jl can be expressed in integral form: 

/ 1 \ \ r M p x - 1 



\XT//" e / dx = y(M) (A4) 

\ Mi 1 1 J x 

where we define a function y(M). It is assumed here that in case of Mi = the corresponding term involving this 
cell does not occur at all in Au, for example if Mi — and M 2 = 2, Au would be Au = (v\ x + v 2x )/M 2 i.e. involves 
only velocities from cell 2. In other words, we assume 1/Mi to be zero for Mi = 0. 
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From Eq. I|A4|) it follows that y(M) satisfies the following linear differential equation 



dM " M 
For M <C 1 we obtain the following approximative solution: 



dy l-e" M /A , 

" 1 y = (A5) 



3 f 1 

y(M) ~ M - -M 2 + — M 3 - ... (A6) 



4 36 

whereas for large M, M>1 the asymptotic solution is 

Let us now consider (|Au|) where the average over the particle numbers cannot be done independently for every 
cell but involves a double sum: 



1 1 \\ 2M v 2 - A T M n+m „ , M ^, /l M n /4oN 
+ - 6 £ \ln + m + 2e £ V n " = ^ (A8) 

n,m — 1 ro=l 

where the second sum comes from the special treatment if one of the two cells is empty. 

For very large M , M% and Mi are simply replaced by their average value M and one obtains 

9(M) ~ ^ (A9) 

For small M just the first few terms in the sums are taken into account and one finds for M< 1: 

g{M) ~ 2M - (4 - 2\/2) A/ 2 + (8 + y / 22/3 - W2) M a + ... (A10) 

Hence, the behavior of (|Au|) j\JksT is propotional to M^FlpK— M 2 (4 — 2\/2)/v / 27r at small M, reaches a maximum 
around M = 1 and behaves as ~ 1/VttM at very large M. This non-monotonic behavior suggests that functions of 
7 should not occur in the equation of state but should be cancelled by an appropriate collision probability. Without 
the cancellation there might be unphysical phase transitions at densities around M = 1. 



[1] T. Ihle, D.M. Kroll, Phys. Rev. E 63, 020201 (2001). 

[2] M. Ripoll, K. Mussawisade, R.G. Winkler, G. Gompper, Europhys. Lett. 68, 106 (2004). 
[3] A. Malevanets, R. Kapral, J. Chem. Phys. 110, 8605 (1999); 112, 7269 (2000). 
[4] T. Ihle, E. Tuzel, D.M. Kroll, Phys. Rev. E 70, 035701 (2004). 

[5] N. Kikuchi, CM. Pooley, J.F. Ryder, J.M. Yeomans, J. Chem. Phys. 119, 6388 (2003). 

[6] T. Ihle, D.M. Kroll, Phys. Rev. E 67, 066705 (2003); 67 066706 (2003). 

[7] M. Hecht, J. Harting, T. Ihle, H.J. Herrmann, Phys. Rev. E 72, 011408 (2005). 

[8] J. T.Padding, A. A. Louis, Phys. Rev. Lett. 93, 220601 (2004). 

[9] Y. Inoue, Y. Chen, H. Ohashi, J. Stat. Phys. 107, 85 (2002). 
[10] A. Malevanets, J.M. Yeomans, Europhys. Lett. 52, 231 (2000). 
[11] H. Noguchi and G. Gompper, Phys. Rev. Lett. 93, 258102 (2004). 
[12] K. Tucci and R. Kapral, J. Chem. Phys. 120 , 8262 (2004). 
[13] T. Sakai, Y, Chen, H. Ohashi, Phys. Rev. E 65, 031503 (2002). 
[14] T. Ihle, E. Tuzel, D.M. Kroll, Europhys. Lett. 73, 664 (2006). 
[15] E. Tuzel, T. Ihle, D.M. Kroll, Math, and Comput . in Simul. 72, 232 (2006). 
[16] E. Tuzel, T. Ihle, D.M. Kroll, cond-mat/0606628 to appear in Phys. Rev. E. 
[17] F.J. Alexander, A.L. Garcia, B.J. Alder, Phys. Rev. Lett. 74, 5212 (1995). 

[18] G.A. Bird, Molecular Gas Dynamics and the Direct Simulation of Gas Flows, (Clarendon Press, Oxford, 1994). 
[19] A.L. Garcia, W. Wagner, Transp. Theory and Stat. Phys. 31, 579 (2002). 

[20] K. Nabu, in Rarefied Gas Dynamics 15, ed. by V. BofR and C. Cercignani (Teubner, Stuttgart, 1986), pp. 369-383. 
[21] J.M. Montanero, A. Santos, Phys. Rev. E 54, 438 (1996). 



29 



[22] A. Frezzotti, Phys. Fluids 9, 1329 (1997). 

[23] E. Tiizel, M. Strauss, T. Ihle, D.M. Kroll, Phys. Rev. E 66, 036701 (2003). 
[24] CM. Pooley, PhD thesis, Oxford University, 2004. 
[25] E. Tiizel, G. Pan, T. Ihle, D.M. Kroll, in preparation. 

[26] J.M. Thijssen, "Computational Physics", Cambridge University Press, 1999. 

[27] R. Zwanzig, in Lectures in Theoretical Physics, Vol. 3, p. 135 (Wiley, New York, 1961). 

[28] H. Mori, Prog. Theor. Phys. 33, 423 (1965). 

[29] H. Mori, Prog. Theor. Phys. 34, 399 (1965). 

[30] J.W. Dufty and M.H. Ernst, J. Phys. Chem. 93, 7015 (1989). 

[31] M.H. Ernst, E.H. Hauge, and J.M.J, van Leeuwen, Phys. Rev. Lett. 25, 1254 (1970). 
[32] M.H. Ernst, E.H. Hauge, and J.M.J, van Leeuwen, Phys. Rev. A 4, 2055 (1971). 

[33] D. Forster, Hydrodynamic Fluctuations, Broken Symmetry, and Correlation Functions, (W.A. Benjamin, Reading, 1975). 



